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Abstract 

This article characterizes certain small multistationary chemical reaction networks. Specif- 
ically we identify the 'smallest CFSTR atoms of multistationarity', namely those containing 
one non-flow reaction, which may be irreversible or reversible. We will refer to such atoms as 
one-reaction atoms of multistationarity. CFSTR atoms of multistationarity were introduced 
in the recent work of Joshi and Shiu. We recall that a fully open network (alternatively a 
fully open CFSTR) is a chemical reaction network where all chemical species participate in 
the inflow and the outflow; and a CFSTR atom of multistationarity N is a multistationary 
fully open network such that the operation of 'removing reactions' or 'removing species' 
destroys the multistationarity of N. The class of CFSTR atoms of multistationarity de- 
termines the entire class of multistationary fully open networks since possessing a CFSTR 
atom of multistationarity as an 'embedded network' (N is embedded in G if -/V can be ob- 
tained from G by removing a finite subset of reactions and species) was proven recently to 
be a sufficient condition for a fully open network to admit multistationarity. 

We find that there are infinitely many one-reaction atoms, however the set of such atoms 
can be characterized using two types, each type containing two parameters. The first type 
contains one chemical species while the second type contains two chemical species; while both 
types contain one irreversible reaction. We identify both types with the chemical process 
of autocatalysis. Furthermore, we give a complete classification by multistationarity of all 
fully open networks which contain one non-flow (possibly reversible) reaction. Moreover, 
we obtain new sufficient conditions for establishing multistationarity of certain fully open 
networks beyond the one-reaction setting. 

Keywords: chemical reaction networks, mass-action kinetics, multiple steady states, 
deficiency one theorem, deficiency one algorithm, atoms of multistationarity. 

1 Introduction 

Chemical reaction networks are used to model systems that occur in chemical engineering and 
systems biology. The property of existence of multiple positive steady states (also known as 
multistationarity) provides a necessary condition for a biochemical network to act as a switch. 
Therefore it is an important problem to determine for which chemical reaction networks there 
exist positive rate constants for which the network has multiple positive steady states. De- 
termining whether a chemical reaction network admits multiple positive steady states is dif- 
ficult: for instance, in the mass-action kinetics setting, it requires determining existence of 
multiple positive solutions to a system of multivariate polynomials with unknown coefficients. 



Several criteria exist which help rule out multistationarity in chemical reaction networks. Im- 
portant examples of such criteria are Deficiency Zero and Deficiency One Theorems of Feinberg 
[3 [B], the Jacobian criterion and the more general injectivity test of Craciun and Feinberg 
[H [21 [31 HI [HI [13]. Such criteria may be thought of as providing necessary conditions for mul- 
tistationarity since avoiding these conditions is necessary for multiple steady states. On the 
other hand, sufficient conditions for multistationarity are relatively rare. Instances where mul- 
tistationarity can be established include Feinberg's Deficiency One Algorithm [8] and Ellison 
and Feinberg's Advanced Deficiency Algorithm [5]. These criteria have been implemented in 
the Chemical Reaction Network Toolbox, software available online for free download and use 

Within the fully open network setting (a fully open network is a chemical reaction network 
where all chemical species participate in inflow and outflow), recent results by Joshi and Shiu 
[12j give a new approach for establishing multistationarity via 'atoms of multistationarity' 
(see Definition 14. 9p . Possessing an atom of multistationarity as an 'embedded network' (see 
Definition I2.2p is a sufficient condition for multistationarity in fully open networks. Using this 
approach, the problem of classifying multistationary fully open networks may be reduced to two 
relatively simpler problems: 1) determining the atoms of multistationarity, and 2) determining 
whether a network possesses one of the known atoms of multistationarity as an embedded 
network. Here we focus on the first problem, and provide an answer for the smallest networks. 

Example 1.1. Consider the following fully open networks Nl — N3 in species A, B, C, D and 
E. In other words, the flow reactions <^ A, «^ B, <^ C , <^ D, and T± E are part of 
all the three networks. 

Nl: A + B -> A + C 2B -»• A + D A + 2E 3E. 
N2: A + B A + C 2B A + D A + E^2E. 
N3: A + C -> A + B 2B -> A + D A + E^2E. 

Which of these three networks admit multiple positive steady states? Note that N2 differs from 
Nl only in the third reaction and in the stoichiometric coefficients of the species E. Moreover, 
N2 differs from N3 only in the direction of the first reaction. However, we will show that 
only Nl and N3 admit multiple positive steady states (both by virtue of possessing atoms of 
multistationarity) . 

In this work, we characterize the class of 'smallest' atoms of multistationarity, namely 
those containing one non-flow reaction, which may be irreversible or reversible. This is a 
continuation of the work in |12j . where the authors catalog all two-reaction bimolecular atoms 
of multistationarity. Atoms of multistationarity containing one non-flow reaction will be referred 
to as one-reaction atoms of multistationarity. Consider the following general one-reaction fully 



2 



open network consisting of s species all of which are in the inflow and outflow: 

k 

+± Xi , f < i < s 

k 

aiXi + a 2 X 2 + ■■■ + a s X s biXi + 6 2 A 2 + ■■■ + b s X s 

k b 

where at least one of the rate constants k a or kb is assumed to be positive. The ki and li 
are positive rate constants which indicate the rate at which the species Xi flows in and out, 
respectively. The stoichiometric coefficients aj and b{ are assumed to be non-negative integers. 
The main theorem (Theorem 14. If) in this article gives a simple arithmetic relation on the 
stoichiometric coefficients which establishes whether the network is multistationary or not. 

As corollaries to Theorem 14.11 we get two- important results. The first result is Theorem 
14.101 which gives a classification of the entire set of one-reaction atoms of multistationarity. We 
find that the infinitely many one-reaction atoms of multistationarity can be classified into two 
types, each type parametrized by two integers. The first type contains one chemical species 
and the second type contains two chemical species. Furthermore, the non-flow reaction in both 
types of atoms is irreversible. In other words, we find that there are no one-reaction atoms of 
multistationarity with a reversible non-flow reaction. The second result is Theorem 14. Ill which 
is obtained by combining Theorem 14.11 with the 'embedded network theorem' of Joshi and Shiu 
|12| and extends the applicability of Theorem 14.11 beyond the setting of one-reaction networks. 
Theorem 14.111 states that a fully open network with any number of non-flow reactions admits 
multiple steady states if it possesses a one-reaction atom of multistationarity as an embedded 
network. 

We find that the multistationary one-reaction fully open networks including the one-reaction 
atoms of multistationarity are identified with the chemical process of autocatalysis. Specifically 
a one-reaction fully open network is multistationary if and only if the network contains a 
non-flow reaction with a set of species that are autocatalytic {i.e. they appear with a higher 
stoichiometric coefficient in the product complex than in the reactant complex) , and the sum of 
the stoichiometric coefficients in such autocatalytic species in the reactant complex is at least 
two. 

Other authors have previously approached the problem of identifying the smallest chemical 
reaction networks with a certain specified property. Smallest multistationary chemical reaction 
networks with the mass-preserving property have been studied in [T3]. The smallest chemical 
reaction outside the fully open network setting (by number of species, number of reactions, 
and number of terms in the differential equation) was studied in [T7] and the smallest chemical 
reaction network with Hopf bifurcation was studied in [181 fl9] . Other examples of classification 
by multistationarity of small networks include [9J [T5| [T6] . Recently generalized catalytic and 
autocatalytic networks have been studied in [TO] . 

This article is organized as follows. Section 2 provides the background information on chem- 
ical reaction networks including the basic definitions, notation and the Deficiency One Theorem 
of Feinberg. Section 3 provides a review of the Deficiency One Algorithm of Feinberg. In Section 
4, we state and prove our main theorem which gives a characterization of one-reaction fully open 
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networks by multistationarity. The first corollary gives a complete classification of one-reaction 
atoms of multistationarity. As a second corollary we get sufficient conditions for establishing 
multistationarity of larger networks by way of one-reaction atoms of multistationarity. 

2 Chemical reaction network theory 

We begin with a review of the notation and basic definitions related to chemical reaction 
networks. An example of a chemical reaction is the following: 

X 1 + 2X 2 -> A 3 

The Xi are called chemical species, and X\ + 2X 2 and A3 are called chemical complexes. 
Assigning the reactant complex X\ + 2X 2 to the vector y = (1, 2, 0) and the product complex 
A3 to the vector y' = (0,0, 1), we rewrite the reaction as y — » y'. 1,2 and are respectively 
the stoichiometric coefficients of the species X\ , X 2 and A3 in the reactant complex X\ + 2 X 2 . 
An autocatalytic reaction is a chemical reaction in which at least one chemical species which 
appears in both the product and the reactant complex has a higher stoichiometric coefficient in 
the product complex than in the reactant complex, for instance A+B — > 2A+C is autocatalytic. 
Many of the definitions in this paper follow those in Joshi and Shiu [12] ; now we define chemical 
reaction networks. 

Definition 2.1. Let S = {X{\, C = {y}, and 1Z = {y — > y'} denote finite sets of species, 
complexes, and reactions, respectively. The triple {S,C,1Z} is called a chemical reaction network 
if it satisfies the following: 

1. for each complex y € C, there exists a reaction in 7Z for which y is the reactant complex 
or y is the product complex, and 

2. for each species Xi £ S, there exists a complex y £ C that contains A,. 

A subset of the reactions 71' C 1Z defines the subnetwork {S\c\ n C\n.i , TV}, where C\n> 
denotes the set of complexes that appear in the reactions TV, and S\c , denotes the set of 
species that appear in those complexes. We now define the notion of an embedded network, a 
more general notion than a subnetwork. 

Definition 2.2. Let (3 = {S,C,1Z} be a chemical reaction network. 

1. Consider a subset of the species S C S, a subset of the complexes C C C, and a subset of 
the reactions R C 7Z. 

• The restriction of R to S, denoted by R\s, is the set of reactions obtained by taking 
the reactions in R and removing all species not in S from the reactant and product 
complexes. If a trivial reaction (one in which the reactant and product complexes 
are the same) is obtained in this process, then it is removed. Also removed are extra 
copies of repeated reactions. 
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• The restriction of C to R, denoted by C\r, is the set of (reactant and product) 
complexes of the reactions in R. 

• The restriction of S to C , denoted by S\c, is the set of species that are in the 
complexes in C. 

2. The network obtained from <8 by removing a subset of species {Xi} C S is the network 

{s\{Xi}, C\ nlsX{Xi} , Tl\ S \{ Xi }} ■ 

3. The network obtained from (3 by removing a set of reactions {y — > y'} C 7Z is the subnet- 
work 

{ 5 lc|„ XW} , Ch\{y^y>}, 1l\{y^y'}} . 

4- Let <3 = {S,C,7Z} be a chemical reaction network. An embedded network of &, which 
is defined by a subset of the species, S = {X^,Xi 2 , . . . ,Xi k } C S, and a subset of the 
reactions, R = {Rj 1 , Rj 2 , ■ ■ ■ , Rj t } C 7Z, that involve all species of S, is the network 
(S, C\r\ s ,R\s) consisting of the reactions R\s- 

We will use the notation € M' 5 ' for the canonical basis vector whose ith component is 1 
and the other components are 0. Accordingly, ej is the vector corresponding to the complex 
Xi. 

Definition 2.3. 1. A flow reaction contains only one molecule; such a reaction is either an 
inflow reaction — > or an outflow reaction — > 0. A non-flow reaction is any reaction 
that is not a flow reaction. 

2. A chemical reaction network is a continuous-flow stirred-tank reactor (CFSTR) if it con- 
tains all outflow reactions e^ — > (for all Xi G S) and a CFSTR is a fully open network 
if it contains all inflow reactions — > e,/. We note that a fully open network is referred 
to as a fully open CFSTR in ITE/ . 



3. A one-reaction fully open network is a fully open network with either one irreversible 
non-flow reaction or one reversible non-flow reaction. 

Let Xi represent the concentration of the chemical species Xi and let x = (x±, . . . , xigi). 
Let yfc = (yki, ■ ■ ■ , Uk\s\) De the stoichiometric coefficients of the species (Xi, . . . , X\$\) in 
the reactant complex of the k-th reaction. For a vector of positive reaction rate constants 
(ki, K2, ■ ■ ■ , G M^q we will assume that the reaction rate is determined by mass-action 
kinetics, or in other words by 

where by convention 0° = 1. 
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A chemical reaction network endowed with mass-action kinetics generates the following 
system of mass-action ODEs: 



\n\ 

x(t) = *£ Kk x(ty*(y' k -y k ) =: f(x(t)) . (2) 

k=l 

We define the function sign on real numbers by 

' 1 if .x > 
sign(x) = < if x = 
k -l ifx<0 

Definition 2.4. i. T/ie stoichiometric subspace of a network is the vector space spanned by 
the reaction vectors of the network, S := span({y' — y|y — > y' € TZ}). A vector fi € M' 5 ' 
is said to be sign compatible with the stoichiometric subspace if there exists a A € S such 
that sign(pi) = sign(Xi), for 1 < % < \S\. Note that ([2]) implies that a trajectory x(t) that 
begins at a positive vector x(0) = c° € M^q remains in the stoichiometric compatibility 
class, which we denote by 

V := (c° + S) n M§ , (3) 

for all positive time; in other words, this set V is forward-invariant with respect to ([2]). 
Two points in the same stoichiometric compatibility class V are said to be stoichiometri- 
cally compatible. Note that y and y' are stoichiometrically compatible ify — y'&S. A 
subset of all stoichiometrically compatible positive vectors forms a positive stoichiometric 
compatibility class. 

2. A concentration vectors € M^q is a (positive) steady state of the system (|2|) z/f(x) = 0. 

A steady state x is nondegenerate i/Imdf(x) = S. (Here, "df(x) " is the Jacobian matrix 

8fi 

d.l: j 



off at x: the \S\ x \S\-matrix whose (i,j)-th entry is equal to the partial derivative 5^-(x) ). 



In the case of a fully open network, the reaction vector for the i-th inflow reaction is the 
f-th canonical basis vector of M) s \, so the stoichiometric subspace is S = ffi' 5 '. It follows that for 
a fully open network, the unique stoichiometric compatibility class is the nonnegative orthant: 
V = 

A chemical reaction network can be thought of as a graph whose vertex set is C and whose 
edge set is 1Z. The next few definitions address the graph-related structure of a chemical reaction 
network. 

Definition 2.5. 1. The complexes y and y' are adjacent if either y —¥ y' G 1Z or y' — s> y € 
1Z. A linkage class of a network is a connected component of the network. 

2. A complex y is said to communicate with the complex y' if there exists a subset of C, 
{yi, . . . ,y„,zi, . . . ,z m }, such that (y -4 yi ->■ y 2 -)• . . .y„_i -4 y n ->■ y') and (y' 
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Zi — > Z2 — > ■ ■ ■ z n _i — >• z n — > y) are reactions in 1Z. A terminal strong linkage class of the 
network is a set of complexes within a linkage class which communicate with each other 
but do not communicate with any complex outside the set. 

Definition 2.6. The deficiency of a chemical reaction network denoted by 6 is defined to be 
5 = n — I — d where n is the number of distinct complexes in the network, I is the number of 
linkage classes and d is the dimension of the stoichiometric subspace. 

The Deficiency Zero Theorem and the Deficiency One Theorem of Feinberg establish that 
a network with deficiency zero and a certain subclass of networks with deficiency one cannot 
admit multiple mass- action steady states. 

Theorem 2.7 (Deficiency Zero Theorem, Feinberg [7]). Suppose that a chemical reaction net- 
work N has deficiency and that each linkage class of N is a terminal strong linkage class. 
For all positive reaction rate constants, the mass-action ODEs have precisely one steady state 
in each positive stoichiometric compatibility class. 

Theorem 2.8 (Deficiency One Theorem, Feinberg [7]). Consider a chemical reaction network 
endowed with mass action kinetics, and with I linkage classes, each containing just one terminal 
strong linkage class. Suppose that the deficiency of the network is 5, that the deficiencies of 
the individual linkage classes are 5j,j = 1, . . . ,1, and that these numbers satisfy the following 
conditions: 

(i) Sj < 1, 

fa) £j=i^- = 6 - 

Then, for an arbitrary choice of rate constants, the chemical reaction network does not admit 
multiple steady states within a positive stoichiometric compatibility class. 

The Deficiency One Algorithm stated in Section [3] requires the following regularity condition 
on the network. Most networks arising as models of chemical processes satisfy this regularity 
condition. 

Definition 2.9. A network is considered to be regular if it satisfies the following conditions: 

1. The reaction vectors of the network are positively dependent. In other words, there exists 
a set of positive numbers {a yi ^ yj \yi — > yj £ 1Z} such that J2n a yi~^yj(yj ~ yO = 0- 

2. Each linkage class in the network contains just one terminal strong linkage class. 

3. For each pair of adjacent complexes {y^y^} in a terminal strong linkage class of the 
linkage class T, if both reactions y» — > yj and y,- — > y^ are removed, then T disconnects 
into two linkage classes. 
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3 A review of the Deficiency One Algorithm 



Here we review the Deficiency One Algorithm of Feinberg [8]. The Deficiency One Algorithm 
takes as input a regular deficiency one network with two or more linkage classes, each of defi- 
ciency zero, and determines whether the network permits multiple steady states or not. The 
algorithm has two variations depending on whether or not the network contains an irreversible 
reaction. We now describe the two algorithms. In the following we will let s := \S\ denote the 
number of species in the network, r := \1Z\ denote the number of reactions in the network and 
fi = Gui, ...,£t s ) G R s . 



3.1 Deficiency one algorithm for a network that contains irreversible reac- 
tions 

Input: A regular deficiency one network with two or more linkage classes, each of deficiency 
zero and such that there is at least one irreversible reaction in the network. 

Step 1. Determine a set of numbers, {g\, g2, ■ ■ ■ , 9\e\} n °t & U zero, such that the following hold: 

( a ) Ei=i 9iYi = 0. 

(b) The gi corresponding to complexes in each linkage class sum to zero. 

(c) The gi corresponding to complexes in each terminal strong linkage class sum to 
a nonnegative number. 

Step 2. For a pair of adjacent complexes y p and y q in a terminal strong linkage class, remove 
the reaction arrows between the pair. Because of regularity Condition 3, the linkage 
class containing this terminal strong linkage class disconnects into two disjoint com- 
ponents. Sum over the gi associated with the complexes in one of the resulting two 
components of the linkage class. Write y p • fi — y q ■ fi > (respectively = 0, or < 0) 
depending on whether sum over of gi is positive (respectively is zero, or is negative). 
Repeat this step on the original network for every distinct pair of adjacent complexes 
in all terminal strong linkage classes. 

Step 3. Partition the set of reactant complexes in the network into three subsets U, M and L 
as follows: 

(a) All complexes that do not belong to a terminal strong linkage class are placed in 
the subset M. 

(b) All complexes in the same terminal strong linkage class are placed in the same 
subset. 

Step 4. For the partition chosen in step 3, for each pair of distinct complexes {yi,yj} in the 
subset M, write the relation yj • fj, = yj • //. 

Step 5. For the partition chosen in step 3, do the following. 
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(a) For each complex yj in U and each complex in M write y$ • fx > yj ■ fx. 

(b) For each complex y» in M and each complex y^ in L write v,- • fj, > yk ■ fJ>- 

(c) For each complex y^ in U and each complex y& in L write yj • fi > yk ■ fJ-- 

Step 6. For the partition chosen in step 3, do the following. 

(a) For each adjacent pair of complexes in each terminal strong linkage class con- 
tained in U, write the inequality from Step 2. 

(b) For each adjacent pair of complexes in each terminal strong linkage class con- 
tained in L, write the inequality from Step 2 with the inequality sign reversed. 

Step 7. Gather all the relations obtained in Steps 4-6 which results in the inequality system 
for the partition chosen in Step 3. 

Step 8. Determine if there exists a nonzero vector [i which satisfies the inequality system corre- 
sponding to the partition chosen and which is sign compatible with the stoichiometric 
subspace of the network. If there does exist such a vector then the network admits 
multiple positive steady states. 

Step 9. Return to Step 3 and choose a new partition of the reactant complexes. Repeat Steps 
4-8, and if necessary repeat this step. 

Output: If there exists a nonzero vector \i which is sign compatible with the stoichiometric 
subspace of the network and which satisfies the inequality system for any partition chosen in 
Step 3, then the network admits multiple mass action steady states. Otherwise, no matter 
what positive rate constants are chosen, the network does not admit multiple mass action 
steady states. 

3.2 Deficiency one algorithm for reversible networks 

Input: A regular deficiency one network with two or more linkage classes, each of deficiency 
zero and such that all reactions in the network are reversible. 

Carry out the same algorithm as the one for a network which contains irreversible reactions. If 
at the end of the algorithm, multistationarity is not established for the network, then repeat 
the algorithm with the signs of the gi chosen in Step 1 reversed. 

Output: Same as Section 13.11 establishes multistationarity or otherwise of the regular defi- 
ciency one network. 

4 Multistationarity in One-Reaction fully open networks 

Now we are ready to state our main theorem. The theorem stated below provides a complete 
characterization by multistationarity of one-reaction fully open networks. The proof of the 
theorem involves an application of Feinberg's Deficiency One Algorithm [8]. 
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Theorem 4.1. 1. Consider a fully open network endowed with mass action kinetics and 
which contains only one non-flow reaction: 

mXi + a 2 X 2 + ■■■ + a s X s -»• 61X1 + b 2 X 2 + • • • + b s X s , 

where a^bi > 0. Then the fully open network admits multiple positive steady states if and 
only if the following holds: 

Oi>l. (4) 

i: bi>Oi 

Moreover, these multistationary fully open networks admit nondegenerate steady states. 

2. Consider a fully open network endowed with mass action kinetics in which the only non- 
flow reactions consist of a pair of reversible reactions: 

atXx + a 2 X 2 + ■■■ + a s X s ±5 biXi + b 2 X 2 + ■■■ + b s X s , 

where ai,bi > 0. The fully open network admits multiple positive steady states if and only 
if the following holds: 



a,i > 1 or bj > 1 ■ (5) 



E 

i: bj>a,j i: a,j>b 



Moreover, these multistationary fully open networks admit nondegenerate steady states. 



Remark 4.2. Theorem \4- 1\ establishes a relation between autocatalysis and multistationarity in 
one-reaction fully open networks. More precisely, conditions @ or ([5]) may be interpreted as 
follows: a one-reaction fully open network is multistationary if and only if it has an autocatalytic 
reaction and the total molecularity of the autocatalytic species in the reactant complex is at least 
two. 

Before we prove the theorem, we need a few technical lemmas. In the following two lemmas, 
we establish the existence of nondegenerate multiple positive mass-action steady states of certain 
simple but important one-reaction fully open networks. 

Lemma 4.3. Let a 2 > 01 > 1. Consider the following fully open network N containing one 
non-flow reaction: 

§ X ai X A a 2 X (6) 

lx 

Let k* := ( — ] ( — ) . Then the following holds: 

a 2 -a 1 \aij V k x J 



Ifk { 



€ (0, k*) N has two nondegenerate positive mass-action steady states 
= k* N has one positive mass-action steady state 
> k* N has no positive mass-action steady states. 
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Proof. Let x represent the concentration of the species X. The network ([6]) when endowed with 
mass-action kinetics results in the following dynamical system: 

x = f{x) = kx — Ixx 4- k(a 2 — a\)x ai 

Clearly /(0) = kx > and f(x) — > oo as x — > oo. Suppose the minimum of f(x) occurs 

at x*. Since f'(x*) = —lx + ka\{a 2 — ai)(x*) ai_1 = 0, we have x* = ( fc ai (a^_ ai ) ) ^ 1 anc ^ 

/(x*) = kx + lx (j^ — lj sc*- If /(^*) < (= 0, > 0) then there are two (one, zero resp.) 
positive steady states. Solving these inequalities for k gives the desired condition. For the case 
where there are two steady states, the nondegeneracy of the two steady states is clear. □ 

Lemma 4.4. Let b\ > 1 and b 2 > 1. The following fully open network M containing one 
non-flow reaction admits multiple nondegenerate positive mass-action steady states. 

0^1 +± Y X + Y A b x X + b 2 Y (7) 

Furthermore, M has two positive nondegenerate mass-action steady states if and only if the 
parameters satisfy the following inequality: 

( lx + hkx(b 2 - 1) - ky(h - 1))V > 1 (8) 



4k(h - l)l x k x V Iy 

In the case of equality in the above equation, M has one positive mass-action steady state, and 
in the case of the reverse inequality M has no positive mass-action steady states. 

Proof. Let x and y represent the concentrations of the species X and Y respectively. The 
network ([7]) when endowed with mass-action kinetics results in the following dynamical system. 



(9) 



x \ = / fi(x,y) \ = f kx ~ Ixx + k(h - l)xy 
y ) V h( x >v) ) V k Y - Ivy + - ^)xy 

Some straightforward calculation reveals that the zeros of 

f(x,y) := (fi(x, y), f 2 (x, y)) coincide with the zeros of the following system 

g{x) := k x - (lx + j^( k xQ>2 - 1) - k Y (h - 1))) x + k(h - X)^-x 2 = (10) 

So that N has two steady states if and only if the second equation in (|1U|) has two positive 
roots which occurs when g(x*) < where x* is the minimum of g(x). The inequality g(x*) < 
is equivalent to ([8]) . It only remains to show that the set of positive parameters satisfying the 

inequality (jH| is non-empty. To this end, let fey ■= kx (j^E^j , k := 2 (bi-i)l x ' an< ^ ^ x = kx + 1- 

With these choices, the left side of ([8|) is ^f^— - which is greater than 1 for all positive kx- 
The nondegeneracy of the two steady states is clear. 

□ 
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Remark 4.5. The reaction networks studied in the Lemmas \4-3\ and \4-4\ are 'one-reaction 
atoms of multistationarity ' (see Definition \4 ■ The lemmas establish that these one-reaction 
atoms admit multiple positive nondegenerate mass-action steady states and that the one-reaction 
atoms admit at most two such steady states. By the theorem of Joshi and Shiu (Lemma \4- 7\ ), 
it follows that any fully open network which contains one of these atoms of multistationarity as 
an embedded network admits at least two positive nondegenerate mass-action steady states. 

Lemma 4.6. Let {ai, 02, ■ ■ ■ , a s , b\, 62, • • • , b s } be a set of nonnegative integers. Consider the 
following system of inequalities: 

s 

^^a-jyUj > max fj,j > (11) 
i=l —3— 

sign(fii) = sign(bi — a^) (1 < % < s). (12) 
This system has a solution fi* € R s \ {0} if and only ifY2vb >cu ai ^ 1' 

Proof. We first assume that the inequality system has a nonzero solution denoted by fi* . So 
the set {j\fJ,j > 0} = {j\bj > aj} is nonempty and 

i:bi>ai i:bi>a,i i:bi<a,i i:bi=a,i 

s 

= } Oi^u* > max /i* = max fi* > 0. 

j:i<i<s r- b 3> a o 

where in the last line we used (llip twice. This shows that 

> aiU* > max n* (13) 

t:bi>ai 

If Yli-b^^ a i — then there exists at most one i such that bj > or > 0. So that Yli-.b >a a *^«* — 
[i~ which contradicts (fT3|) . So we must have X^-ft^oj a * > ^ 

Conversely, assume that Yli-b >a a » - > 1- ^ or au * sucn that 6j > Oj, choose /ij = 1 and for 
all j such that bj < aj, choose fij = —e for some e > 0. This choice clearly satisfies (|12p and 

s 

^ a-ilii = ^ v-i^i + X] aifjii 

i=l i:bi>a,i i:bi<ai 

= ^2 a i — e ^2 a « — ^ — e a i > 1 = m ^x ^ > 

i:bi>ai i:bi«ii i:bi«ii 

where the last inequality follows by choosing e sufficiently small. □ 

If a fully open network TV" is an embedded network of a fully open network G, then we can 
extend the steady states of N to G using the following result. 
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Lemma 4.7. [Theorem 4-2 and Corollary 4-6 of Joshi and Shiu HE/ J Let N be a fully open 
network embedded in a fully open network G. 

• If N admits multiple nondegenerate positive mass-action steady states, then so does G. 
Moreover, if N admits finitely many such steady states, then G admits at least as many. 

• If N admits multiple positive exponentially stable steady states then so does G. Moreover, 
if N admits finitely many such steady states, then G admits at least as many. 

The following lemma deals with some simple cases of one-reaction fully open networks 
appearing in the statement of Theorem 14. 1\ which can be handled by using either Deficiency 
Zero Theorem or Deficiency One Theorem of Feinberg [7] . 

Lemma 4.8. Consider the following one-reaction fully open network N: 

0<^ej 1 < i < s 
y a ^ Yb 



where y b — > y a may have a zero reaction rate 
is possibly irreversible). If {y a ,yb} H {0,ei, . . . 
action steady states. 



constant (in other words, the non-flow reaction 
, e s } / then N does not admit multiple mass- 



Proof. We will show that if {y a , y&} n {0, ei , . . . , e s } ^ holds then the network has a deficiency 
of either zero or one. So that either by applying Deficiency Zero or Deficiency One Theorem, 
multiple mass-action steady states can be ruled out. Since we are assuming that y a — > y b is a 
non-flow reaction, it is not the case that one of the complexes {y a >yfc} is the complex and 
the other is a unimolecular complex. So suppose that one of the complexes (either y a or y b but 
not both) in the non-flow reaction is either the complex or is unimolecular. Then the network 
has only one linkage class (1 = 1) and n = s + 2 complexes so that the deficiency of the network 
is 5 = n — I — d = (s + 2) — (1) — (s) = 1. If the reaction y b — > y a has a positive rate constant, 
then the network is reversible and the entire unique linkage class is a terminal strong linkage 
class. Otherwise, the unique linkage class has exactly one terminal strong linkage class, either 
y b or C \ y a , depending on whether y a G {0, ei, . . . , e s } or y b € {0, ei, . . . , e s }, respectively. In 
either case, the hypotheses of the Deficiency One Theorem (Theorem I2.8|) are satisfied which 
rules out multiple steady states for such a network. 

On the other hand if both y a and y^ are unimolecular, then iV has one linkage class (I = 1), 
and n = s + 1 complexes. So the deficiency of N is S = n — I — d = (s + 1) — (1) — (s) = 0. 
Whether the rate constant for the reaction y^ — > y a is zero or positive, the entire unique linkage 
class is a terminal strong linkage class. So by the Deficiency Zero Theorem (Theorem I2.TH . N 
does not admit multiple mass-action steady states. □ 

Now we are ready to prove the one-reaction fully open network theorem. 

Proof of Theorem \4-l\ Assume first that {y a , y&}n{0, ei, . . . , e s } ^ 0. This means that at least 
one of the complexes y a or y b is at most unimolecular, which implies that 6 >a . ai < 1 and 
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Si- a,i>bi — 1- Furthermore by Lemma [4.81 the one-reaction network does not admit multiple 
positive mass-action steady states. 

So it suffices to assume that {y a ,yb} H {0,ei, . . . , e s } = 0. From hereon, we will assume 
it to be the case that each of y a and y^ is at least bimolecular. In the following we will let 
C a := a\X\ + a 2 X 2 + • • • + a s X s and C b '■= b\X\ + 62-^2 + • • • + b s X s and the corresponding 
complex vectors are y a = (a±, . . . , a s ) and yj, = (61, . . . , b s ). For the general one-reaction fully 
open network, 

a. The set of species is S := {Xi\l <i < s}. 

b. The set of complexes is C = {e^l < i < s} U {0} U {Ylt=i Yll=i bi&i}- We relabel the 
complexes y» = ej for 1 < i < s, y s+ i = 0, y s+2 = y a = (01, • • • ,a s ) and y s+3 = y b = 
(61, . . .,b s ). 

c. The set of reactions is either 1Zi rrev = {0<^ej|l < i < s} U {y a — > y;,} or 1Z rev = {0 i=i 
e»|l < i < s} U {y a — > y;,} U {y^, — > y a } depending on whether the non-flow reaction is 
irreversible or reversible. The two connected components of the network partition the set of 
complexes into C\ := {ej|l < i < s} U {0} and C2 ■= {J2t=i a i&i-> Ylt=i bi&i}- 

The proof will proceed by application of the Deficiency One Algorithm. In 1 through 4 in the 
following, we lay the groundwork for application of the Deficiency One Algorithm by checking 
the conditions of validity of the algorithm. In 5, we apply the Deficiency One Algorithm to 
the one-reaction fully open network containing an irreversible non-flow reaction. In 6, we apply 
the Deficiency One Algorithm to a one-reaction fully open network containing a reversible non- 
flow reaction, and finally in 7, we establish the nondegeneracy of the steady states for the 
multistationary fully open networks. 

1. Dimension of the stoichiometric subspace. For a fully open network, the stoichio- 
metric subspace has 'full' dimension. For s species, the dimension of the stoichiometric 
subspace is d = s. 

2. Deficiency of the network. Since we are assuming that neither of the two complexes 
in the non-flow reaction is the complex or is unimolecular, {y a ,yb} H {0, ei, . . . , e s } = 0. 
The number of complexes is n = s + 3, the number of linkage classes is I = 2 and the total 
deficiency is 8 = n — I — d = (s + 3) — (2) — (s) = 1. On the other, the deficiency of the 
two linkage classes is given by 

(a) {0 ^ Xi\l < i < s}: <5i = m - 1 - d x = (s + 1) - (1) - (s) = 0. 

(b) For {y a f± y 6 }: 8 2 = n 2 - 1 - d 2 = (2) - (1) - (1) = 0. 

Such networks satisfy the hypotheses of the Deficiency One Algorithm [8] provided they 
are also regular. In the following we check the regularity of the network. 

3. Regularity of the network. 
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(a) Regularity Condition 1: Note that the outflow and inflow reactions for species j have 
the reaction vectors e,- and —ej. The set {ej, € S} forms a basis over the non- 
negative integers, which implies that the reaction vectors are positively dependent. 

(b) Regularity Condition 2: The entire linkage class {0,ei, . . . ,e s } is a terminal strong 
linkage class. Either the linkage class {y a ,yb} 1S a terminal strong linkage class or 
contains the terminal strong linkage class {y&} depending on whether the non-flow 
reaction is reversible or not. 

(c) Regularity Condition 3: If the inflow and outflow reactions of some species j are 
removed or if the non-flow reaction were to be removed, the linkage class containing 
that reaction will be disconnected into two linkage classes. 

This shows that the network under consideration satisfies the hypothesis of the Deficiency 
One Algorithm. 

4. Sign compatibility with the stoichiometric subspace. Since the network under 
consideration is a fully open network, the stoichiometric subspace has "full" dimension. 
In other words, the stoichiometric subspace S = W. Every vector fx* = (//i, . . . , /i s ) 
is contained in S and therefore is trivially sign compatible with S. Thus it suffices to 
determine an inequality system for a partition obtained in Step 3 of the algorithm and 
then to determine if the inequality system has any solution. If there is such a solution for 
any partition, then the network admits multiple steady states. 

5. Application of the Deficiency One Algorithm to a network containing at least 
one irreversible reaction. Suppose first that the non-flow reaction is not reversible. 
In other words (y& — > y a ) ^ 1Z. The set of reaction vectors for the network is 1Zi rrev = 
{ej, — &i\i € S} U {Yli=i{bi — flj) ej}. Let r := cardinality (1Z) = 2s + 1. 



Step 1. Let 



(1 < i < s) 

s 

9s+l ■= ~ 9s+2 ■■= -1, 9s+3 ■■= 1- (14) 



1=1 



We now check that the gi satisfy the three required conditions. 

\C\ s / s \ s s 

(a) ^9iYi = ^2(ai-bi)ei+{ - a;) + (-1) ^ a& + (1) ^ b& = 0. 

i=l i=l \i=l / i=l i=l 

(b) It is clear that the gi sum to zero for the two linkage classes C\ = {e^ 1 1 < 
i < s} U {0} and C 2 = {£- =1 a&i, £- =1 b t e,}. 

(c) C\ is a terminal strong linkage class while {^i=i ls a terminal strong 
linkage class within The corresponding gi sum to and 1, respectively. 

Step 2. C 2 does not contain a pair of complexes in its terminal strong linkage class. In 
C±, all adjacent pairs of complexes are of the type {0, e^} for some i. For a fixed 
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k, let y p := {e^} and y q := 0. So the sum over gi in y p is gt = &k — bk and 
y p ■ fi — y q ■ fi = efc • \i — = fik . So that for all k we get the following system of 
inequalities: 

sign(/ii) = sign(oj - fcj) (1 < i < s) (15) 

Step 3. The set of reactant complexes is C \ {y&}. 

(a) y a is the only complex that does not belong to a terminal strong linkage 
class and so y a € M . 

(b) C\ is a terminal strong linkage class and so all the complexes in C\ must be 
placed in the same subset. The three choices of where to place the complexes 
in C\ lead to the following three partitions of the reactant complexes. 

(i) 17 = 0, M = {yJU{0}U{e i |l<i<s}, 7 = 0. 

(ii) 17 = 0, M = {y a }, L = {0}U{e l \l<i<s}. 

(iii) U = {0} U {e 4 jl < i < s}, M = {y a } 7 = 0. 

It is straightforward to show that if the inequality system resulting from partition 
(ii) has a solution fi* , then the inequality resulting from partition (iii) has a 
solution which is equal to — /x*. This is true whenever two partitions are related 
by switching the contents of the sets U and 7 [8] . Thus we will restrict attention 
to partitions (i) and (ii). 

Step 4. For partition (i), we get the following system of equations. 



H = (1 < i < s), ^2 cum = 0, 

i=i 

s 

^2 aiPi = V-j (1 < 3 < s), Iii = lij (1 <i,j < s) 



i=l 

which clearly does not have a nonzero solution (j,* = (//*, . . . ,/i*). So it suffices 
to consider partition (ii) only. Beginning here and in all the following steps, we 
will assume that we are considering partition (ii) even when this is not explicitly 
stated. For partition (ii), there is only one complex in M and so we do not get 
any equations from applying this step. 

Step 5. Since there is one complex in M and s + 1 complexes in 7, we get the following 
system of s + 1 inequalities 



Oi/it > o, ^2 a '^i > n (i < i < s ) 



i=\ i=l 



Step 6. For each adjacent pair of complexes in 7, we write the inequality from Step 2 
with the inequality reversed. 

sign(/ij) = sign(6j - a/) (1 < i < s) 
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Step 7. We gather all the inequalities from Steps 4-6 to get the inequality system. 



s s 

1=1 t=l 

sign(/ij) = sign(6j - m) (1 <i< s). 

The first inequality holds only if there exists a j such that fij > 0. So the first 
s + 1 inequalities are equivalent to X)i=i > maxi<j< s > 0. The system 
of inequalities to be solved may be written as: 

s 

y diHi > max fij > 
i=i 

sign(^) = sign(6j - m), (1 < i < s) (16) 
Step 8. By Lemma 14,61 the system of inequalities (|16p has a solution if and only if 

Hi:bi>ai a i > 1- 

Step 9. We have already shown in Steps 3 and 4 that we only need to consider partition 
(ii) since the other two partitions (i) and (iii) do not provide any new information. 

This completes the proof of the theorem for networks that contain at least one irreversible 
reaction. 

6. Application of the Deficiency One Algorithm to a reversible network. Now 

suppose that the non-flow reaction is reversible which results in a reversible network, in 
other words a network containing only reversible reactions. The set of reaction vectors 
for the network is TZ rev = {e i: - e,|i € S] U {X)i=i(&i ~ U {X^i=i( a i ~~ h)^i}- 

The three networks TZ rev = {0 +± e^l < i < s} U {y a -> y b } U {y b -> y a }, lZ irrev = {0 z± 
e»|l < i < s} U {y a ->■ y b }, and 7^i rrei , = {0 |1 < i < s} U {y^ ->■ y a } generate the 
same stoichiometric subspace. By Lemma 14. 7| if either 1Zi rrev or lZi rre v admits multiple 
positive nondegenerate mass-action steady states then so does 1Z rev . This shows that if 
either £V 6 >a . a« > 1 or Xh-<j>&- ^* > ^> then 7^ rei , admits multiple steady states. 

Assume now that Y2i-b >a- °* — ^ an< ^ Si-a >6 ^ — 1- It only remains to show that the 
network does not permit multiple steady states. Note that for reversible networks we 
need to go through the steps of the algorithm twice - once for the gi chosen according to 
Step 1, and then again for gi which are negative of those chosen according to Step 1. 

Step 1. For 1 < i < s + 2, we let gi be the same as in (|14p . Since the complexes and 
the linkage classes are the same as in the irreversible case, there is nothing to 
check in parts (a) and (b). For part (c), we note that for a reversible network 
every terminal strong linkage class is a linkage class and so the gi corresponding 
to complexes in all terminal strong linkage classes sum to zero. 
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Step 2. We obtain the same set of inequalities as in (|15p for the terminal strong linkage 
class C\. Further we note that 1 = g s +z > g s +2 = — 1, so we get the following 
system of inequalities: 

sign(jLti) = sign(aj - h) , 1 < i < s 

s 

22(k - Oi)fM > 
i=i 

Step 3. The only constraint is that all complexes in C\ belong to the same subset and all 
complexes in C2 belong to the same subset. Thus there are 9 distinct partitions. 
However, Remark 4.1.G in Feinberg [8] tells us that we only need to examine 
partitions for which both subsets U and L are nonempty, since the condition that 
all terminal strong linkage classes contain more than one complex is satisfied. 
Furthermore, as in the irreversible case, interchanging contents of the subsets U 
and L result in equivalent inequality systems in the sense that if one inequality 
system has a solution then so does the other. Thus we only need to consider 
the following partition U = C 2 = "Ei=i a i&u Yli=i M = 0, L = C\ = 

{0} U {ej|l < % < s}. 

Step 4. Since M = 0, no equations result from this step. 

Step 5. Comparing complexes in U and L we get the following system 



^ a^i > 0, ^2 a '^i > N 0- < 3 < s ) 

i=l i=l 

s s 

kfJ-i > 0, ^2 > jjbj (1 < j < s). 



i=l i=l 

Step 6. We write the inequality from Step 2 with the inequality corresponding to the 
complex in L reversed. 

sign(jUi) = sign(6j - ca) , 1 < i < s 
s 

i=l 

Step 7. We gather all the inequalities from Steps 4 to 6. 

s s 

bifn > ^2 a i^i > max Mj > 

1=1 1=1 —j— 

sign(/ii) = sign(6j - Oj), (1 < i < s) (17) 
Step 8. By Lemma f4.6l this system of inequalities has no nonzero solutions for Y^i-b >a ai — 
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As mentioned earlier, for the reversible case, we need to carry out the algorithm again 
with signs of all gi reversed. This results in an interchange of the roles of a% and bi and 
so it is straightforward to see that we get the following system of inequalities in Step 7. 

s s 

ciiHi > bifJ-i > max Mj > o 

i=l i=l —3 — 

sign(/Xj) = sign(aj - 6;), (1 < i < s) (18) 

Since by hypothesis Y2i- a >b- — 1> it follows from Lemma 14,61 that the above system of 
inequalities does not have a nonzero solution. This completes the proof of the theorem 
for the reversible case. 

7. Nondegeneracy of the steady states for the multistationary fully open net- 
works. We have shown that every multistationary one-reaction fully open network has 
an embedded fully open network of the form satisfying the hypotheses of either Lemma 
14.31 or of Lemma 14.41 Since the fully open networks appearing in Lemma 14.31 and Lemma 
14.41 admit multiple positive nondegenerate mass-action steady states, by Lemma 14.71 it 
follows that if one of the one-reaction fully open networks admits multiple positive mass- 
action steady states, then it admits multiple positive nondegenerate mass- action steady 
states. 

This completes the proof of the theorem. □ 

We recall the following definition of CFSTR atom of multistationarity of Joshi and Shiu 
EE]. 

Definition 4.9. 1. A fully open network is a CFSTR atom of multistationarity if it admits 
multiple nondegenerate positive mass-action steady states and it is minimal with respect 
to the embedded network relation among all such fully open networks. 

2. A one-reaction atom of multistationarity is a CFSTR atom of multistationarity containing 
one non-flow (irreversible or reversible) reaction. 

3. A fully open network G is said to possess a CFSTR atom of multistationarity if there 
exists an embedded network N of G that is a CFSTR atom of multistationarity. 

As a corollary of Theorem 14.11 we get a complete classification of all one-reaction atoms of 
multistationarity. We state this result as a theorem. 

Theorem 4.10. 1. A one-reaction fully open network is a CFSTR atom of multistationarity 
if and only if it consists of one irreversible non-flow reaction and that non-flow reaction 
has one of the following two forms: 

axX -> a 2 X , or X + Y -> h X + b 2 Y , (19) 

where a% > a± > 1, or, respectively, b\ > 1 and b 2 > 1. 
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2. A one-reaction fully open network possesses a CFSTR atom of multistationarity as an 
embedded network if and only if it admits multiple nondegenerate positive mass-action 
steady states. 

Proof. Evidently, the one-reaction fully open networks in (I19p admit multiple nondegenerate 
positive mass-action steady states by Lemmas 14.31 and 14.41 Furthermore, it is clear that both 
types of fully open networks are minimal in the class of multistationary fully open networks 
with respect to the embedded network relation. 

On the other hand, assume that N is a one-reaction fully open network with multiple steady 
states. Then by Theorem 14. 11 N has a subnetwork containing a reaction of the following form: 

ciXi + c 2 X 2 + ■■■ + c s X s ->• di.Xi + d 2 X 2 + ■■■ + d s X s 

with Ylvdi>a c i > 1- If there exists a j such that dj > Cj > 1, then we let X = Xj, a\ = Cj, 
and a 2 = dj. Otherwise, there exists a pair of indices (i,j) such that Cj = Cj = 1, di > 1 and 
dj > 1. In that case, we let (X,Y) = {X i: Xj), b\ = di, and b 2 = dj. 

The second part of the theorem follows from Lemma 14.71 This completes the proof. □ 

The one-reaction atoms of multistationarity are useful beyond the one-reaction setting. If a 
one-reaction atom of multistationarity is an embedded network of a fully open network N with 
possibly more than one non-flow reaction, then the 'embedding theorem' of Joshi and Shiu 
(Lemma 14. 7p . may be used to infer that N has multiple positive nondegenerate mass-action 
steady states. We state as a theorem the following important corollary of Theorem 14.11 whose 
scope of application is beyond the one-reaction setting. 

Theorem 4.11. A fully open network (which may contain more than one non-flow reaction) 
admits multiple nondegenerate positive mass-action steady states if it contains as an embedded 
network a reaction of one of the following forms: 

1. a±X —?■ a 2 X with a 2 > ai > 1. 

2. X + Y -> b\X + b 2 Y with bi>l and b 2 > 1. 

Recall from |12j that bimolecular networks are such that each complex in the network has 
at most two molecules. In other words, a complex in a bimolecular network has one of the 
following forms: 0, A, 2 A or A + B. Theorem |4J] establishes that there are no bimolecular one- 
reaction fully open networks with multiple steady states. This in turn implies that the smallest 
(by number of reactions) bimolecular fully open networks with multiple steady states should 
contain at least two non-flow reactions (where one reaction is not merely the reverse of the other 
reaction). In fact, there do exist bimolecular two-reaction CFSTR atoms of multistationarity 
and this set has been fully catalogued in [T2] , 

We end by posing the question of identifying and characterizing 'larger' CFSTR atoms of 
multistationarity, i.e. the ones that contain more than one non-flow reaction. Since tests that 
involve checking whether a certain large atom of multistationarity is embedded in a network 
may be computationally difficult, characterization of the CFSTR atoms of multistationarity in 
terms of general principles may be particularly helpful. 
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Example 4.12 (Example 11.11 continued) . We end by answering the question posed in Example 



1. For network Nl the fully open network Gl containing the non-flow reaction 2E —¥ 3E 
is an embedded network. Since Gl is a one-reaction atom of multistationarity, it follows 
that Nl is multistationary. 

2. For the fully open network N3, first remove the reaction A+E —¥ 2E, and then remove the 
species C and D . This gives the fully open network G3 containing the non-flow reactions: 
A — > A + B and 2B — > A. G3 is one of the known two-reaction bimolecular atoms of 
multistationarity U2\j and therefore N3 is multistationary. 

3. A little bit of calculation shows that N2 possesses no known atoms of multistationarity. 
In fact, plugging the fully open network N2 into the Chemical Reaction Network Toolbox 

reveals that N2 does not admit multiple positive steady states. 
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